
A single grid local mode analysis is used to predict the smoothing properties of numerical schemes 
for solving the Navier-Stokes equations with factorization based on Stone’s Strongly Implicit Method. 
Four difference approximations for the convection terms are considered, namely, hybrid, central, second- 
order upwind, and third-order upwind. Smoothing factors from the analysis are compared with practical 
convergence factors in a multigrid method for flow over a backward facing step and it is found that the 
local mode analysis correctly predicts the effects of Reynolds number and higher-order schemes. 

1 INTRODUCTION 


The successful use of multigrid methods to accelerate convergence rates is dependent on the ability of 
the numerical algorithm to dampen high frequency error components since these components cannot be 
resolved on coarser grids. High frequency components have short coupling ranges; therefore, their 
smoothing is a localized process meaning that only one isolated computational stencil need be analyzed 
and the effect of boundaries can be neglected. This is the approach of local mode analysis for the predic- 
tion of smoothing properties which was first introduced by Brandt [1] for various partial differential equa- 
tions and numerical algorithms. Shaw and Sivaloganathan [2] extended this analysis to the SIMPLE 
pressure correction algorithm using alternating direction implicit (ADI) relaxation for the solution of the 
algebraic system of equations for varying Reynolds numbers and under-relaxation factors. Convection 
terms were approximated using a hybrid of first-order upwind and second-order central differencing. 

The present paper uses local mode analysis to predict the smoothing properties of numerical algo- 
rithms for calculation of two-dimensional recirculating flows using higher-order difference schemes for 
convection terms introduced via deferred correction and Stone’s Strongly Implicit Method for factoriza- 
tion of the resulting system of algebraic equations. Reynolds number and higher-order convection 
approximation effects are addressed and compared to multigrid results for laminar flow over a backward 
facing step. 


PRECEDING PAGE BLANK NOT FILMED 



2 THEORETICAL ANALYSIS 


2.1 Governing equations 


The equations governing steady, two-dimensional, incompressible flow can be written as: 


«L<p») + £(pv) = 

3(p,v) + A(pv J ) --g+P 



where u and v are the velocity components in the x and y directions, respectively, p is the pressure, q 
is the absolute viscosity, and p is the density. Equations (1) - (3) represent the conservation of mass and 
momentum in the x and y directions, respectively. 


The solution sequence is a predictor-corrector method which follows the SIMPLE algorithm of Patan- 
kar and Spalding [3]. Factorization of the system of equations is based on Stone’s Strongly Implicit 
Method. The flow geometry and boundary conditions are shown in Figure 1. 

The governing equations are discretized by integrating over a set of three staggered control volumes 
and the locations of variables are shown in Figure 2. The central control volume is for the pressure. Equa- 
tion (2) can be discretized by integrating over the left-shifted control volume for the u component of 
velocity. This leads to: 
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Equation (3) is discretized by integrating over the bottom-shifted control volume for the v component 
of velocity: 
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where the a, coefficients contain convection and diffusion terms, the subscripts of u, v, and p refer to 
the location of the variables (see figure 2), h x and hy are the grid spacing in the x and y directions, respec- 
tively, and p o is the absolute viscosity of the fluid, assumed constant 


2.2 Approximation of convection terms 

The di coefficients in equations (4) and (5) are dependent on the approximation used for the convec- 


tion terms. The present analysis investigates the hybrid, central, second-order upwind (2nd OU), and 
third-order upwind (3rd OU) approximation schemes. 


The coefficients for the hybrid scheme have the following form: 

"jt 

V a |P u oL r pM ° m 

a w = maxiO,-^-^} + max{-^, 0} 
n x x 

f r, |P V ol A f ^ V ° fU 

ntax{0, -t-tjt-} +max\ , 

h) 2h y n y 

max{ 0.^-^p) +max{j^,0} 
h) 2h y y 

a EE ~ a ww = °NN = a SS = 0 
a P = £a, 


a N = 

o s 


a P - 

is taken over the a coefficients, u 0 and v 0 are the frozei 
' equations (2) and (3), the max{a,b} operator selects the 
esents the absolute value. 

ing general form: 


( 6 ) 

(7) 

( 8 ) 

(9) 

( 10 ) 

(ID 


where the sum for cip is taken over the a coeffi' 
due to the linearization of equations (2) and (3), the 
ments a and b, and 1 1 represents the absolute value. 

The coefficients for the higher-order schemes have the following 
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The values of// and/ 2 depend on the higher-order method to be used. For central differencing, f } =f 2 
= 0, for 2nd OU differencing,// = 1/2 and/? = 1 . and for 3rd OU differencing,// = 1/8 and/? = 1/4. 


2.3 Local mode analysis 


The Strongly Implicit Method (SIP) by Stone [4] solves the algebraic equations shown in equations 
(4) and (5) in a fully implicit manner. All variables are treated as unknowns, as opposed to line-relaxation 
methods which consider only lines of constant x or y as unknowns while sweeping through the computa- 
tional domain. In the local mode analysis outlined below, the variables that are updated at the end of a 
relaxation sweep will be denoted by a dot over the variable, such as u P , while those from the beginning of 
the relaxation sweep or those unchanged by the current relaxation sweep will be written as above, such as 
Up. Higher-order approximations for the convection terms are introduced via deferred correction (see 
Khosla and Rubin [5]). In this procedure, the a t coefficients are calculated initially using equations (6) - 
(11) for the hybrid scheme. As the solution proceeds, the higher-order scheme is slowly introduced via 
corrections to the source terms. At the end when the solution is fully converged, the coefficients are effec- 
tively those of the higher-order scheme outlined in equations (12) - (20). The base hybrid coefficients will 
be denoted by an additional subscript h , such as a P l A , while the higher-order coefficients will not have an 
additional subscript and will be written as above, such as a p . 

The relaxation of equation (4) (the discretized x momentum equation), with under-relaxation and 
deferred correction can be written as: 
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where r m is the under-relaxation factor for the x and y momentum equations, ry is the relaxation factor 
for introducing higher-order coefficients and is set to unity for the analysis. The exact solution for U y V y 
and P also satisfies equation (21). If an equation written with the exact solution for U, V y and P is sub- 
tracted from equation (21), which is written in terms of the approximate solution, u, v, and p , the error in 
the solution can be introduced. The error has components defined as; e u = U-u y e v = V-v, and 
& = P-p. 

Equation (21) written in terms of the error becomes: 
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Since the continuous governing equations (1) - (3) have been linearized during the discretization, a 
single Fourier component of the error can be considered as. 
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where i = 0, and 0 2 are the components of the phase angle vector, a e , which is the error ampli- 

tude of the single Fourier mode ©,, 0 2 . Similar expressions exist for other grid points to the east and 
north and for the variables v and p. Substituting the single Fourier modes into equation (22) and dividi g 
through by ^ [e > x ^ h * +0 2 > y M 9 equation (22) becomes: 
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where: s, = sin (0^2) 

s 2 = sin (0 2 /2) 


Equation (24) can be written in a more compact form, if the following variables are defined: 
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Equation (24) can then be written as: 


a« = l(va“ 9 -<pa v e-CV e ) 


(25) 


Following the same procedure for equation (5) (the discretized y momentum equation) yields: 
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where: 
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Equations (25) and (26) can be combined to give the amplification matrix A j for the complete opera- 
tion. This yields: 



In compact form: 


«e = ( 27 ) 

Equation (27) yields the amplification matrix defining how the amplitude of the Fourier mode, with 
phase angles 0j and 0 2 , is amplified during relaxation of the x and y momentum equations. 

The SIMPLE pressure correction of Patankar and Spalding [3] follows the relaxation of the x and y 
momentum equations. A single dot over a variable will denote a value at the completion of the relaxation 
of the x and y momentum equations. A double dot over a variable will denote a value at the completion of 
the pressure correction. The variables u, v, and p are corrected following Shaw and Sivaloganathan |2|: 


U P~ u ( 5 Pp 5 Pw) 
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where r uv is the relaxation factor for correcting u and v velocities and r„ is the relaxation factor for 
updating pressure. The value 8 p is a pressure increment such that the velocity field u and v will satisfy 
conservation of mass. It is obtained by discretizing equation (3) and substituting for the velocities and 
corrections given in equations (28) - (30). This yields an equation for the pressure correction: 

cfp§p P = d > N&P N + d > s&P s + c f E?>PE + d’w§P w - j-(u E -u P ) - ^-{v N -Vp) ( 31 ) 
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Equations (28) - (31) can be written in compact form as: 
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If equation (35) is solved for bp h and used in equation (32)-(34): 
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This assumes that the pressure correction equation has been solved exactly. As before, the error is 
introduced by writing equations (36)-(38) using the exact solution and then subtracting the result from 
equations (36)-(38) respectively. The errors become: 




Op 

[u 

a P 




i\ = £'’* + r p P \ (8 X h i U h + & h t h ) 


(39) 

(40) 

( 41 ) 


The Fourier components of the error can be substituted as before to give the A 2 amplification matrix 
which governs the amplification of errors during the pressure correction phase: 
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In compact form: 


“a = 4 2 6c e (42) 

where: P h = a F E e‘ 9> + ^ w e~' 9 ' + <f N e 82 + - cf P 

Equation (42) gives the amplification matrix defining how the amplitude of the Fourier mode with 
phase angles 0j and 0 2 is amplified during the pressure correction phase of the algorithm. The amplifica- 
tion matrix for relaxation of the x and y momentum equations and the pressure correction is obtained by 
combining equation (27) and (42) as: 


a e - A 2 A l a % - 4a 0 (43) 

2.4 Smoothing factor 


The smoothing factor is a measure of the worst reduction of the high frequency error components for 
one complete relaxation sweep. It is calculated as the largest eigenvalue of the amplification matrix A 
given by equation (43) for the Fourier modes 0j and 0 2 in the high frequency range defined as: 
ti/2 < |0,| < n and n/2 < 1 0 2 | < n. 


3 RESULTS 


To test the predictive capability of the LMA presented in Section 2, flow over a backward facing step 
(BFS) was computed using a multigrid code based on the FAS-FMG (full approximation storage - full 
multi-grid) algorithm proposed by Brandt [1]. Higher-order schemes were introduced through deferred 
correction only on the finest of three grids with constant grid spacing in the x and y directions. The grid 
sizes from coarsest to finest grid are, n x xn y = 66 x 18, 130 x 34, and 258 x 66, where n x is the number of 
grid points in the x direction and n y is the number of grid points in the y direction. Smoothing properties 
on the two coarser grids are identical for the four schemes since hybrid coefficients were used on these 
grids. Local mode analysis was used to estimate the smoothing factor on the finest grid and this result was 
compared to the number of work units to reach convergence for the multigrid result. The work units 
(WU) and convergence factor (CF) are indicators of the smoothing properties of the algorithm and 
numerical scheme. The work units for a two-dimensional problem with grid refinement in the x and y 
directions are defined as: 

WU = (44) 

1 = 1 

where t is the number of iterations on the i lh grid at convergence, i = 1 for the coarsest grid, and 1 = 
N for the finest grid. The convergence factor is defined as: 

CF = (r/r,) (45) 

where r ( is the initial norm of the residuals of the x momentum, y momentum, and pressure correction 
equation on the fine grid, /y is the norm of the residuals at convergence on the fine grid, and A WU is the 
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change in work units on the fine grid. 


The BFS flow was solved for Reynolds numbers (based on the upstream channel height) of 100, 250, 
and 400 using the hybrid, central, 2nd OU, and 3rd OU schemes. The work units for convergence and 
convergence factors are shown in Table 1 . The smoothing factors were calculated for conditions identical 
to the finest grid of the multigrid results (the same grid spacing, relaxation factors, density, etc..). The 
only parameter varied was the frozen velocity components needed to calculate the coefficients in equa- 
tions (6) - (20) of the various schemes. The maximum velocity components provide an upper bound for 
the smoothing factor which will dominate the smoothing properties since it was found that as the cell- 
Reynolds number (Reynolds number with the length scale based on the grid spacing) approaches zero, 
corresponding to regions where the velocity approaches zero, the smoothing factor also decreases. An 
estimate of the maximum velocity components for the BFS flow is u Q = 1.5 at y = 0.25 at the inlet, and v Q 
= 0.15 near recirculation regions. Three principal flow directions are considered with velocity compo- 
nents given by: u 0 , v 0 = (0, 0. 15), (1.5, 0.15), and (1.5, 0). The SIP method exhibits symmetry about the x 
and y axis so that other flow direction results can be obtained from the three principal flow direction 
results. For example, the smoothing factors for u 0 , v 0 = (-1.5, 0.15), (-1.5, -0.15), and (1.5, -0.15) are 
equal to the smoothing factor for u 0 , v 0 = (1.5, 0.15). The smoothing factor is then defined as the largest 
eigenvalue of the amplification matrix A, defined by equation (43), for the three flow directions while 
restricting the phase angles to the high frequency range. Results from the three principal flow directions 
show that the flow direction u 0 , v 0 = (1.5, 0.15) produced the largest eigenvalue for all Reynolds num- 
bers, and thus the smoothing factor was based on this flow direction. The computed smoothing factors are 
shown in Table 2. 

The results of Table 2 show that as the cell-Reynolds number in the LMA is increased, the smoothing 
factor also increases for the four schemes. More work units will be required to smooth the high frequency 
error components. The results in Table 1 show that the work units increase and the convergence factor 
deteriorates as the Reynolds number increases. For the Re = 100 results, the LMA predicts that the 
smoothing properties of the hybrid, central, and 3rd OU will be virtually identical while that of the 2nd 
OU will be slightly worse. The multigrid results confirm this prediction. For the Re = 250 and 400 results, 
the LMA predicts that the hybrid difference scheme will have the best smoothing properties while the 
central difference scheme will have the worst, and the 2nd OU and 3rd OU difference schemes should be 
similar with the 3rd OU difference scheme slightly better. The multigrid results confirm these predictions 
with the exception being that the 2nd OU difference scheme results converged in slightly less number of 
work units when compared to the 3rd OU difference scheme. Their convergence factors are similar. 


Table I: Work Units/Convergence Factors of Multigrid Results 


Difference 

Scheme 

Re = 100 

Re = 250 

Re = 400 

Hybrid 

59/0.868 

137/0.930 

321/0.981 

Central 

58/0.873 

166/0.964 

569/0.993 

2nd OU 

63/0.891 

148/0.952 

421/0.988 

3rd OU 

58/0.875 

152/0.956 

428/0.988 
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Table II: Smoothing Factors from Local Mode Analysis 


Difference 

Scheme 

Re= 100 

Re = 250 

Re = 400 

Hybrid 

0.902 

0.924 

0.931 

Central 

0.910 

0.950 

0.968 

2nd OU 

0.931 

0.946 

0.950 

3rd OU 

0.899 

0.926 

0.935 

Re &x/R e \ y 

11.72/0.47 

29.30/1.18 

46.88/1.88 


CONCLUSION 


Local mode analysis was performed using four schemes for the approximation of convection terms: 
hybrid, central, second-order upwind, and third-order upwind, over a range of cell-Reynolds numbers. 
The smoothing factors from this analysis were compared with actual multigrid results for flow over a 
backward facing step to test the predictive capability oFIocal mode analysis. It was found that this analy- 
sis is useful in predicting the smoothing properties of the four schemes along with the effect of flow Rey- 
nolds number. This analysis could be extended to predict optimum relaxation factors, grid aspect ratios, 
and other solution algorithms. 
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